The dynamics of simple gene network motifs snbject to extrinsic flnctuations 
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Cellular processes do not follow deterministic rules; even in identical environments genetically 
identical cells can make random choices leading to different phenotypes. This randomness origi¬ 
nates from fluctuations present in the biomolecular interaction networks. Most previous work has 
been focused on the intrinsic noise (IN) of these networks. Yet, especially for high-copy-number 
biomolecules, extrinsic or environmental noise (EN) has been experimentally shown to dominate 
the variation. Here we develop an analytical formalism that allows for calculation of the effect of 
EN on gene expression motifs. We introduce a new method for modeling bounded EN as an auxil¬ 
iary species in the master equation. The method is fully generic and is not limited to systems with 
small EN magnitudes. We focus our study on motifs that can be viewed as the building blocks of 
genetic switches: a non-regulated gene, a self-inhibiting gene, and a self-promoting gene. The role 
of the EN properties (magnitude, correlation time, and distribution) on the statistics of interest 
are systematically investigated, and the effect of fluctuations in different reaction rates is compared. 
Due to its analytical nature, our formalism can be used to quantify the effect of EN on the dynamics 
of biochemical networks and can also be used to improve the interpretation of data from single-cell 
gene expression experiments. 

PACS numbers: 87.16.Yc, 02.50.Ey, 05.40.-a, 87.17.Aa 
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I. INTRODUCTION 

Biochemical processes in cells are inherently noisy, be¬ 
cause many molecular species such as genes, RNAs, and 
proteins that make up intracellular reaction networks 
are present in low copy numbers inside a cell, see e.g. 
Refs. [Il[2]. One of the primary insights to emerge from 
studies on stochastic gene expression is the distinction 
between intrinsic noise (IN) and extrinsic noise (EN) [3]- 
[H]. Experimentally, EN is quantified using the correlation 
in fluctuations between two copies of an identical reporter 
gene expressed separately in the same cell. IN arises from 
fluctuations that are independent for each reporter [^. 
Within a cell, then, IN is the variance due to the dis¬ 
creteness of biomolecules and the probabilistic nature of 
chemical reactions, while EN is the variance arising from 
the fact the genes share a common environment, the cell. 

Such noise in cellular reactions can have important 
consequences, e.g., on cellular decision making. In fact, 
noise can drive cells between distinct gene expression 
states corresponding to different decision phenotypes nni- 
ng. The ultimate stability of a decision state is then 
determined by fluctuations of mRNA and proteins, as 
well as other cellular components, during gene expres¬ 
sion [linHn]- These fluctuations can give rise to spon¬ 
taneous switching between the states, with a switching 
time that depends on their strength and the switch’s ar¬ 
chitecture. Genetic switches can regulate diverse deci¬ 
sion making processes such as microbial environmental 
adaptation, developmental pathways, nutrient homeosta¬ 
sis, and bacteriophage lysogeny [IHHII]- 
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In recent years there have been numerous theoretical 
studies on genetic switches driven by IN, noise arising 
from within a closed system of interest [22lI35j . Yet, 
these and most other studies of genetic switching, have 
neglected sources of EN due to interactions with other 
components in the cell and the environment. EN does 
not arise from a single well-defined process, but rather 
results from the complex chain of events that gave rise to 
a particular cellular state. Variation in a cell’s number 
of ribosomes, transcription factors, and polymerases or 
fluctuations in the cell division time, as well as environ¬ 
mental fluctuations, can all affect the rates of a genetic 
process. These fluctuations in the reaction rates may 
dramatically affect the protein’s statistics including its 
mean, variance and copy-number distribution. Impor¬ 
tantly, in living cells a comparison of the relative con¬ 
tribution of EN versus IN to the protein distributions 
width has shown that EN dominates above copy num¬ 
bers of 0(10 — 100) [53H55] . 

Theoretically, EN has been shown to induce bistabil¬ 
ity [HI I35M3] , vary the distribution tails [5] , and modify 
switching times [43] . In signaling, EN limits the infor¬ 
mation transduction capacity of the pathways [H] 03]. 
It has been shown that EN is present at multiple time 
scales during protein production in bacteria SUIT] and 
that negative feedback can filter EN [15]. Previously, in 
Ref. [48| the authors have studied the interplay between 
IN and EN noise in a genetic switch near bifurcation us¬ 
ing a Fokker-Planck approximation and showed that EN 
can dramatically affect switching. However, their method 
could not be directly used to study questions regarding 
population heterogeneity in metastable systems nor with 
non-Gaussian EN statistics. 

Genetic switches and other more complex circuits using 
mnltiple positive and negative feedback links form the ba- 
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sis of much of the transcriptional regulatory logic in bac¬ 
teria. The overall regulatory network of these microor¬ 
ganisms is commonly composed of repeated patterns of 
relatively small circuits, called motifs. For a review, see 
. As a prerequisite to understanding the effect of EN 
on complex regulatory networks, we would first like to 
understand the role of EN on simple genetic motifs. 

In this work we develop an analytical formalism that 
allows for the quantification of the effect of EN on 
intrinsic-noise-driven gene expression circuits. We intro¬ 
duce a new method for modeling bounded EN as an aux¬ 
iliary species in a master equation that fluctuates accord¬ 
ing to non-Gaussian statistics. We then analyze three ge¬ 
netic motifs: a non-regulated gene, a self-inhibiting gene, 
and a self-promoting gene. These three motifs represent 
the simplest possible circuits, yet commonly occur in bac¬ 
terial transcription networks. We study the properties of 
each motif as related to EN. All our analytical findings 
are tested and compared against numerical Monte Carlo 
simulations. 


II. THEORY 


A. Gene expression under intrinsic noise 


Our starting point is a simple gene-expression model 
without extrinsic noise (EN). Let us denote by n the 
protein’s copy number and by Y 3> 1 the typical pro¬ 
tein abundance in the steady state. Our gene-expression 
model will consist of two reactions: production of pro¬ 
teins at a rate of F{n) and degradation with rate vn. 
Here we assume that the mRNA lifetime is short com¬ 
pared to the cell cycle and momentarily ignore the mRNA 
fluctuations, which will be accounted for in the following. 

In the deterministic picture, the rate equation for the 
protein concentration x = n/N reads 

x = f{x)-x, (1) 

where f{x) = F(n)/N and we have rescaled time by 
the protein degradation rate v. To account for intrin¬ 
sic fluctuations due to the probabilistic reactions and the 
discreteness of the proteins, we write down the chemi¬ 
cal master equation for P„(t) - the probability to find n 
proteins at time t 

Pn = F{n - l)Pn-i + {n + l)P„+i - [F{n) + n]P„. (2) 


We look for the stationary PDF such that = 0. This 
yields a set of recursive equations, whose solution can be 
found analytically. The solution reads na 


TO -I- 1 


m—0 


El- 

m—O 


F{m) 
m + \ 


(3) 


where Pq is a normalization factor such that “ 

1 . 


The stationary solution of Eq. ([^ can also be found 
by using a dissipative WKB approximation [5I1I52]. To 
this end, we assume n ^ 1, treat n as a continuous 
variable, and search for P„ as P„ = P[x) ~ ^-ns{x)_ 
Here, A^ 1 is assumed to be a large parameter, and 
S{x) is called the action. Plugging this ansatz into 
the stationary master equation [Eq. ([^ with P„ = 0], 
we arrive in the leading 0{N) order at a stationary 
Hamilton-Jacobi equation H[x,S'(x)] = 0 with Hamil¬ 
tonian H{x,px) = f{x){e^^ — 1) + x{e~^^ — 1), where we 
have denoted the associated momentum by = S'{x). 
While the trivial zero-energy trajectory Px{x) = 0 of this 
Hamiltonian corresponds to the deterministic dynamics, 
in the leading order the statistics of interest are encoded 
in the nontrivial zero-energy trajectory of this Hamilto¬ 
nian [^, which reads in this case 

Px{x) =\n[x/f{x)\. (4) 

This allows us to calculate the action by integration 
S{x) = J px{x’)dx'. Thus, the PDF and its variance 
due to IN, = NS"{x^)~^ [53], are found to be: 




Nx* 


I - /'(x*)’ 
(5) 


where a:* is the steady-state solution of Eq. Q , and the 
normalization was done over the Gaussian part of the 
PDF around a:*. By transforming the sum into an in¬ 
tegral, Eq. (§ coincides in the leading order in >> 1 
with the PDF in Eq. ([^. 


B. Gene expression under intrinsic and extrinsic 
noise 

Next, we add to our model EN, which is commonly de¬ 
fined as intercellular variability due to fluctuations during 
gene expression that equally affect all genes within a cell. 
We thus introduce EN in the form of one or more fluc¬ 
tuating parameters. For concreteness we assume, e.g., 
that cell-to-cell variability in transcription and transla¬ 
tion rates causes the protein degradation rate v to fluctu¬ 
ate so that V —)• v{t) = ^(t). (In Appendix G we consider 
other fluctuating parameters as well.) As a result, the 
degradation rate becomes n^[t) where ^(t) is a stochas¬ 
tic variable satisfying (^(t)) = 1. Many measured protein 
distributions appear to be well-fit by a negative binomial 
(or gamma) distribution. Without experimental knowl¬ 
edge of how rates fluctuate in vivo, we simply take ^(t) to 
have a negative binomial statistics, as if being controlled 
by a single protein. In addition, ^(t) has variance a^x 
correlation time Tc, satisfying 

Other statistics are also possible in Appendix D 

we consider Ornstein-Uhlenbeck noise. Note that our 
choice of negative binomial statistics for the EN ensures 
that the rates are always positive. 

To model EN, we need a circuit that generates an aux¬ 
iliary species whose copy number fluctuates with nega- 
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tive binomial statistics and correlation time Tc- To cre¬ 
ate one, we use an auxiliary mRNA-protein circuit where 
mRNAs are transcribed at a rate o/tc and degrade with 
rate w/tc, while proteins are translated at a rate uj^/tc 
and degrade at a rate I/tc, which ensures that the cor¬ 
relation time of the auxiliary proteins is Tc- As a result, 
the master equation describing the probability to find m 
auxiliary mRNAs and k auxiliary proteins satisfies: 

Pm,k — (-fm—l,fc Pm.k) “b [(^ “b 

Tc Tc 

H- {Pm,k-1~ Pm,k) H-[(^ + ^)Pm,k+l~ kPm,k]- (6) 

'^C '^c 

As shown in Appendix A, in the limit of short-lived 
mRNA such that w 3> 1, the stationary PDF of the aux¬ 
iliary protein is [54] : 


Tja + k) ( 13 Y ( 1 

r(A: + l)r(a) 'v/3 + i; U + 1 


where Pk is the probability to find k auxiliary pro¬ 
teins. Here k = K^, where K = af3 is the PDF 
mean, while the variance is ^"(1 -b /3). Therefore, choos¬ 
ing a = Ijicr^x ~ 1/-^) P = Kal^ — 1 we find 
that (^) = aP/K = 1 and the variance of ^ becomes 
Ar(l -b P)/K'^ = CTex s-s required by our EN stochastic 
variable. Note, that in the limit of large K such that 
P = — 1 ~ Kal^ and a ~ the negative bino¬ 

mial distribution can be well approximated by a gamma 
distribution Pk ~ P~°‘/T{a) [55]. 

To study the interplay between IN and EN, we combine 
the EN dynamics [Eq. §] with the underlying intrinsic 
noise dynamics [Eq. ([^]. This leads to a 3D master equa¬ 
tion describing the evolution of the probability Pn,m,k to 
find n proteins, m auxiliary mRNAs and k auxiliary pro¬ 
teins, where the death rate of the protein of interest de¬ 
pends on the auxiliary protein. To this end, by using the 
WKB theory and by adiabatically eliminating the short¬ 
lived auxiliary mRNA degree of freedom (see Appendix 
A), we arrive at a stationary Hamilton-Jacobi equation 
p{ = 0 with a reduced Hamiltonian for the protein of 
interest and auxiliary protein: 


= f(x)(eP^ - 1) + —(e - 1) 
P 


+ - 


1 


Ptc [l + Pil-ePi) 


- 1 


+ ( 8 ) 

Tc 


Here we have defined a rescaled EN variable ^ 
where p = K/N is the abundances ratio of the auxil¬ 
iary protein and protein of interest. Also, Px and p^ are 
the momenta associated with the protein of interest and 
the auxiliary protein with corresponding concentrations 
X = n/N and ^ = k/N, while a and P are defined above. 
Hamiltonian ([^ encodes the stochastic dynamics of a 
protein when its degradation rate fluctuates due to nega¬ 
tive binomial EN generated by another auxiliary protein. 
Note, that while the (arbitrary) copy number K of the 


auxiliary protein enters Hamiltonian Q, it does not en¬ 
ter the results below for the statistics of the protein of 
interest. 

Hamiltonian Q can be theoretically analyzed by writ¬ 
ing down the corresponding Hamilton equations, see 
Eqs. (HI) in Appendix B. These can be solved numer¬ 
ically for arbitrary Tc, see Methods section. Analytical 
progress can be made in two important limits: short- 
correlated “white” EN, and long-correlated “adiabatic” 
EN. 


In the white-noise limit, Tc ^ 1, one arrives at a re¬ 
duced white-noise Hamiltonian, which effectively takes 
into account the short-correlated EN. Solving the cor¬ 
responding Hamilton-Jacobi equation we find (see Ap¬ 
pendix B) 


Px = In 


2f{x) V 


1 — VtcX+^ {VtcX— 1)'^ + AV f{x)i 


(9) 


Here, V = is the ratio between the relative EN 

and IN variances (where the IN variance is taken in the 
non-regulated case). Erom Eq. we can calculate the 
action S{x) = Px{x')dx', while the PDF is given by 
Eq. ([^. The action function S{x) cannot be explicitly 
calculated without specifying f{x). Yet, a general re¬ 
sult for the PDF variance can be derived. Differentiat¬ 
ing px{x) [Eq. ([^] once and plugging x = x^ such that 
/(cc*) = a;*, we find the observed PDF variance 


fj 


2 _ 
ohs 


NS"{x,)-^ 


A^a;*(l -I- x^^Vtc) 
1 - /'(a;*) 


( 10 ) 


Comparing with Eq. ([^, this indicates that short- 
correlated EN increases the variance by a factor of 
1 -I- x*Frc. 

In the adiabatic limit, Tc ^ 1, we can assume that the 
EN is almost stationary [IS]. As a result, the protein 
PDE can be written as 


/ OO 

P{0P{n\m, (11) 

-OO 


where P{n\^) is the conditional probability to find n pro¬ 
teins given noise magnitude and P(^) is the probability 
to find EN magnitude In Eq. 0 we effectively opti¬ 
mize the “cost” of reaching a state with n proteins given 
EN magnitude ^ against the probability of choosing such 

e m- 


For simplicity, here we take gamma-distributed EN, 
PiO = ,d““/r(a) where a ~ l/Cg,,, and 

P = P/K ~ This distribution has a mean of 1 and 
variance as required, and is a good approximation 
to the negative binomial distribution for large K [55] . 
Performing the integration in (11) via the saddle-point 
approximation, see Appendix B, the PDF in the adia¬ 
batic limit reads 


P{x)^ 


C 




\/%$[x,^ = C*(a;)] 


( 12 ) 


















4 


where 

In-—dy+ ---, (13 

4(c) /(y) ^ 

is the cost function. Here, ^*(a;) is the solution of the 
saddle point equation 9^<i)(a;,^) = 0, cc = g{^) is the 
dependent stable fixed point found by solving the equa¬ 
tion f{x) = x^, and C is a normalization factor such that 
N J^P{x)dx = 1. 

Similarly as in the white-noise case, here we can also 
calculate the variance of the PDF explicitly for any pro¬ 
duction rate f{x). After some algebra, see Appendix B, 
we find the observed variance of the PDF 


Nxi, 


^obs 


1 - 


1-b 


Vx^ 


1 - /'(a;*) 


(14) 


This indicates that adiabatic EN increases the vari¬ 
ance compared to the IN-only case ([^ by a factor of 
1 + - f'{x*)]. 

Eqs. (10) and (fldl) for the variance are among our main 


results here. When EN is put in the production rate in¬ 
stead of the degradation rate the results for the variance 
in both the white- and adiabatic-EN cases remain the 
same, see Appendix C. 


III. RESULTS 

A. Unregulated gene expression 



FIG. 1. (Color online) Comparison of theory (lines) and 
numerics (symbols) for the non-regulated gene model with 
N = 100 and = 0.2. (a-l-b) Probability distributions for 
white (a; Tc = 0.1) and adiabatic (b; Tc = 1000) EN. Dotted 
lines show the Poisson distribution for the model with only in¬ 
trinsic noise, (c) Observed variance vs EN strength for white 
(red) and adiabatic (blue) noise, (d) Observed variance vs 
EN correlation time for = 0.2. The left curve shows the 
white noise theory and the right shows the adiabatic theory. 


We begin with a model for protein transcription given 
a constant birth rate, namely an unregulated gene. Here 
the rate equation is given by Eq. Q with /(x) = 1, while 
the protein PDE is = e~^ 

Now, we add EN to the protein’s degrada¬ 
tion as described above. In the white-noise 
limit, Tc <C 1, integrating over the momen¬ 
tum ([^ with /(x) = 1, the action function be¬ 
comes S{x) = 1/(Utc) {D(x) -I- ln[UrcX — 1 -I- D(x)] 

+VtcX [In ((x/2)(l — VtcX + D(x))) — 1]}, with 

D(x) = 4 {VtcX — 1)^ -I- 4Utc. Using S'(x), we find 
the PDE, given by Eq. (§, see Figure 0a). Inter¬ 
estingly, in the presence of EN the far right tail of 
the PDF behaves as a power law. Indeed, taking the 
X ^ 1 limit of the PDF, we find a power-law depen¬ 
dence in the leading order P(x) ~ {2VtcX)~^^‘'^^‘=\ 
in contrast to an exponential tail in the IN-only case. 
Nevertheless, because the power-law behavior appears 
only at X >> X* = I, the corresponding probabilities 
are vanishingly small and thus, observing this behavior 
experimentally or even numerically is impractical. 


The variance of the PDF due to white EN [Eq. (10)] 
becomes = N{1 + Vtc) (Figure 0c, d)). That is, 
white EN increases the width of the PDE by a factor of 
VI + Utc. 

In the adiabatic limit, Tc ^ 1, the PDE is given by 
Eq. ( [T^ (Figure0b)). Here, the cost function [Eq. (13)] 


becomes 4'(x,,f) = {i - In^ - 1)/U -b j^f^\n{yC)dy, 
where x = g(^) = 1/^ is the solution to the equation 
/(x) = ^x with /(x) = 1. In addition, the saddle point 

is found at = (1 — Ux)/2 ^1 -b ^1 -b 4U/(1 — Ux)^, 
while (9{4 <I>(x, ^) = (2U -b ^ — Ux^)/(U^^). 


Plugging /(x) = 1 and x* = 1 into Eq. (14), the vari¬ 
ance due to adiabatic EN is given by = N{1 + V) 
(Figure 0c, d)). That is, adiabatic EN increases the 
width of the PDE by a factor of vTbAU which can be 
significant when V > 0(1)- 

To test our theory we performed Monte Carlo sim¬ 
ulations of the full master equation describing all three 
species: the protein of interest n, the auxiliary mRNA m, 
and the auxiliary protein k, see Methods. Figurej^shows 
example comparisons for N = 100 and a ex ranging up 
to 0.5, which are typical values obtained from single-cell 
Escherichia coli protein distributions [36) . Good agree¬ 
ment is obtained between our theory and stochastic sim¬ 
ulations for both the white and adiabatic cases, even for 
quite strong EN. We have also verified that the results 
hold for EN in the birth rate, see Appendix C and Eigure 
S3 in [^. Interestingly, when EN arises in the degra¬ 
dation term, the PDF mean shifts to the right due to 
the nonlinear dependence of the fixed point on the death 
rate. While it is negligible for weak and moderate EN, 
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this shift in the mean becomes significant for very strong 
EN, when aex = 0{\), see the end of Appendix B for 
details. In turn, this shift affects the IN of the system, 
see Appendix E. In such a case IN and EN cannot be 
independently separated, as is commonly assumed (see 
Figure S3 and Figure S4 in [56]). 


B. mRNA-protein model with no feedback 


Now we consider the more realistic case of an unregu¬ 
lated gene but with mRNA present in the model. Here 
mRNAs are transcribed at a rate a, decay with a rate 
7 , and translation of proteins occurs with rate ■jb while 
degradation of proteins occurs with rate 1. As in the aux¬ 
iliary circuit, we take the ratio between the mRNA and 
protein degradation rates to be large 7 ^ 1. The mean 
protein number here is N = ab. The rate equations de¬ 
scribing the average mRNA and protein concentrations, 
r = l/N and x = n/N (with I and n being the respective 
copy-numbers of mRNA and proteins), are f = a/N — "fr, 
and X = Iryr — x. 

In the limit of short-lived mRNA, 7 ^ 1 , the stochas¬ 
tic dynamics has been analyzed by various authors m- 
Using the WKB approximation one can find the protein 
PDF, see Appendix A, which coincides with the negative 
binomial distribution in the limit of n 3> 1. In particu¬ 
lar, the PDF variance becomes N{1 + b), indicating that 
mRNA noise increases the variance by a factor of 1 -I- 6 
compared to the protein-only case [Ti] . 

We now proceed to calculate the observed variance of 
the proteins of interest under negative binomial adiabatic 
EN in the protein’s degradation rate. We do so along 
the same lines done for the protein-only case. Here, ac¬ 
counting for mRNA noise, the momentum given noise 
magnitude ^ becomes Px{x,^) = ln[(I -|- b)x^/{l + bx^)], 
which reduces to the protein-only case for & —>■ 0. In¬ 
tegrating over the momentum, we find the action to be 
S{x,^) = a; In [(I -|- b)x^/{l + bx^)] — 1/(6^) ln(l -|- bx^). 
Now, similarly as done in Eq. ( |T^ , we can define the cost 
function $(x, ^ = S{x, 0 - -S”!^), C] + (C - InC - 1)/^- 
Therefore, the variance of the PDF can be found using 
Eq. ( B1I[ ), by repeating the calculations along the same 
lines as done in Appendix B for the protein-only case. As 
a result, we find the observed variance of the proteins of 
interest, while accounting for mRNA noise, to be 


^obs — N{1 + b + V). 


(15) 


each parameter set. We fit the PDFs to a gamma dis¬ 
tribution to obtain estimates of the gene expression pa¬ 
rameters a fit and bfn- To calculate the accuracy with 
which a fit and bfit recovered the actual parameters, we 
calculated the relative error as Err{a) = \a — afit\/a and 
Err{b) = \b — bfit\/b using the known a and b values from 
the simulations. As can be seen in Figure [^a-|-b), using 
a gamma distribution resulted in poor estimates even in 
the case of relatively weak EN. Given the sensitivity of 
the error to EN, gamma distribution estimates of gene 
expression parameters should be used with caution. 

We then instead used Eq. (151 along with the gene- 
by-gene EN values of to estimate a fit and bfu val¬ 
ues from the simulated dataset. With the mean of the 
distribution N = ab and the observed variance (Jobs = 
N{l + b + V) = ab^l + b + abcr"^^), one can solve directly 
for a = /[(Jlbs-N-N'^al^) and b = 

Using the calculated mean and variance of the PDF, 
along with the known from the simulations, we re¬ 
covered estimates for a fit and bfu- Figure ^a-|-b) shows 
that if one knows the strength of the EN for a gene, 
Eq. (15) can reliably recover the true a and b values until 
the total variance becomes dominated by EN for V > 10. 

Next we attempted to recover the a fit and bfit values 
from a genome-scale protein abundance data set from E. 
coli [36] . Here, we made a simplifying assumption that 
a constant global EN of Gex = 0.31 influenced all genes 
equally (see Figure S5 in |5Bj). We estimated the afit 
and bfit values for each gene using the gamma distribu¬ 
tion and also using Eq. (15) with this global EN. Figure 
[^c-f) shows a comparison of the two methods. A few 
trends are apparent from the results. When accounting 
for EN, the global saturation in the burst frequency a dis¬ 
appears and instead we see a continuous linear increase 
in the burst frequency. Likewise, an observed global in¬ 
crease in b values at higher V disappears and a more 
uniform distribution of b values is seen with respect to 
V. Since V is correlated with overall expression levels 
(IN goes down as V goes up) this implies that burst fre¬ 
quency is a significant driver of protein expression levels 
in E. coli. Figure S6, see [5S], shows the fits versus mean 
expression. Gene-by-gene estimates of EN, rather than a 
single global EN, would lead to even better estimates of 
gene expression parameters. 


C. Self-inhibiting gene 


The gamma distribution is widely used to analyze 
single-cell protein abundance data |SS]. The protein’s 
PDF is fit to a gamma distribution and the a and b val¬ 
ues resulting from the fit are interpreted as the gene’s 
burst frequency and burst size, respectively. We wanted 
to study how EN would affect such interpretations. To 
this end, we performed a large number of stochastic sim¬ 
ulations across a wide range of values for a, b, and a^x, 
again for biological ranges seen in single-cell experiments, 
and calculated stationary PDFs using 10^ data points for 


Next we consider the case of a self-inhibiting gene. A 
self inhibiting gene is a simple yet common motif that is 
capable of filtering some types of IN m, although it can 
lose effectiveness when multiple time scales are involved 
[5H] . We were therefore interested to study the ability of 
a self-inhibiting gene to filter EN. 

The rate equation is given by Eq. (IH where we took 
the production rate to be f{x) = {1 + ^ /{1 + fix), and f 
is the inhibition strength. Here we chose a simple form 
of nonlinear inhibitory Hill-like function with Hill coef- 
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X Fit w/ gamma 
oFit w/ Eq. (15) 




10' 

V 



bfitw/ Eq. (15) 


the PDF is given by Eq. ([^. Using Eq. ( |Io| ) with 
f{x) = (l+/3)/(l+^a;) anda;* = 1, the observed variance 
in this case becomes 


— 

^obs 


N{l+P) 
1 + 2/3 


(l + ^Te). 


(16) 


This result indicates that negative inhibition can elimi¬ 
nate EN. Indeed, returns to its non-regulated value 
without EN, TV, when the inhibition strength satisfies 
/3 = VtcI{1 — Vtc), which holds as long as Vtc < 1. 
That is, our choice of negative inhibition with h = 1 
can only attenuate moderate EN, and can reduce the ob¬ 
served variance at most by a factor of 2. 

In the adiabatic limit, we can find the PDF using 
Eqs. ((I2| and (fl^ (see Appendix B for details), with 
g{0 = l/(2/?0(-?+Ve"+ 4/3(/3 + l)e). Using Eq. (M 
with f{x) = (1 -F /?)/(! -F /3x) and /'(s* = 1) = 
—f3/{f3 + 1), the observed variance is 


^obs 


Nil +13) 
1-F2/3 


' , vq + py 

1 + 2/3 _ 


(17) 


Figure l^a-c) shows good agreement between theory and 
simulations over a wide range of parameters. Again, we 
see that negative inhibition can eliminate EN when /3 = 
(Vl + 4U -F 2U — l)/[2(2 — U)]. Here, the maximum 
EN that can be attenuated for this particular choice of 
inhibition is U = = 2. 


D. Higher order inhibition 


FIG. 2. (Color online) (a-Fb) Relative error in afu and hja, 
respectively, from simulations spanning a wide range of pa¬ 
rameters. Blue crosses show the values from htting to a 
gamma distribution. Red circles show hts from Eq. ( |15| .(c-Fd) 
The a Jit and bju values vs V from htting the genome-scale 
protein abundance data from using (blue crosses) the 
gamma distribution or (red circles) Eq. (151 with a global EN 
of (Tea: = 0.31. (c-Ff) Relationship between aju and bjit val¬ 
ues obtained from a gamma distribution and Eq. (151. Points 
are colored by logio{V) where V = The solid y = x 

lines are a guide to the eye. 


ficient h = 1 (below we will consider higher values of 
h as well), whose fixed point a;* = 1 coincides with the 
non-regulated gene. 

To find the PDF in the IN-only case, we integrate over 
Eq. (|^ using/(x) = {1 +(3)/{1 +/3x). This yields 5'(x) = 
—2x -F (1//3) ln(l -F f3x) -F xln [x(l -F /3x)/(l -F /?)], while 
the PDF is given by Eq. ([^ with x* = 1. The PDF 
variance NS"{x^)~^ = N{1 + /3)/{1 + 2/3), indicates that 
such negative inhibition decreases the PDF variance by 
a factor of 2 at most. 

Adding negative binomial EN into the protein degra¬ 
dation rate, in the white noise limit, the momentum is 
given by Eq. ([^ with /(x) = (1 -F /3)/(l + I3x), while 


We now consider a more generic inhibition function 
/(x) = (1 -F /3)/(l + /3x^) with arbitrary Hill-coefficient 
h. Here, in the white-noise limit, we find 


^obs 


N{l+/3) 

1 -F /3{h + 1) 


{1 + Vt,). 


(18) 


In this case, EN can be eliminated by taking /3 = 
VTc/{h — Vtc), which holds for Vtc < h. In the adia¬ 
batic limit we obtain 


N{l + /3) 

1 -F /3{h + 1) 


' , v{i + P) 

^ l + /3{h + l) 


(19) 


Here, EN is eliminated when /3 = {hy'l + AV + 2V — 
h)/[2{h{h -F 1) — V)], which can be achieved as long as 
V < Vmax = h{h+l). One can see that as h is increased, 
this inhibition mechanism becomes more efficient in elim¬ 
inating EN. 

We wanted to examine the relationship between the 
critical inhibition strength /3cr that will exactly elimi¬ 
nate EN and the inhibition order. We calculated /3cr 
across h values ranging from 0.1 to 100 for various values 
of V. Our analytical framework provides a significant 
advantage over simulations for studying such large pa¬ 
rameters spaces. Figurej^d) shows that cooperativity in 
inhibition is a necessary feature for systems that dampen 
strong EN. 
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FIG. 3. (Color online) Comparison of theory (lines) and 
stochastic simulations (symbols) for the self-inhibited gene 
model with N = 100 and adiabatic EN Tc = 100. (a) Proba¬ 
bility distributions for (solid line) theory and (x’s) numerics 
with j3 = Per = 1.09 and h — 3. Dotted line shows the Pois¬ 
son distribution for the model with only intrinsic noise and 
triangles show the distribution in the absence of negative in¬ 
hibition. (b) Observed variance vs EN strength with /3 = 1.0 
and various values of h. (c) Observed variance vs inhibition 
strength P for aex = 0.2 and for the same h values as in (b). 
(d) The critical P and h values that exactly cancel EN for a 
given relative EN strength V = 


E. Self-promoting gene 

Lastly, we consider the case of a self-promoting gene. 
Self-promoting genes can serve as genetic switches, allow¬ 
ing the cell to change between two alternate expression 
states. The rate equation now satisfies Eq. 0 where 
we take the production rate to be /(a;) = ao + (l — 
ao)9{x — xq) — a step function imitating a Hill-like func¬ 
tion with a high Hill coefficient of positive feedback [15] . 
Here ao < xq < 1, where ao is the protein concentra¬ 
tion in the off state, xq is the threshold concentration, 
and N is the protein abundance in the on state such that 
X = n/N. Unlike previously [48] . our derivation here is 
generic and does not require xq to be near one of the 
metastable states ao or 1. 

In the absence of EN, the mean switching time (MST) 
from the off to the on states and vice versa, Toff^on 
and Ton^off, can be calculated using the master equa¬ 
tion Q and employing the WKB approximation (see e.g., 
Ref. P8| ). Indeed, assuming that we start from the vicin¬ 
ity of the off or on metastable states we find the mo¬ 


menta Poff(x) = ln(a:/ao) and Pon{x) = ln(a;), where 
we have used the fact that f[x) = ao for x < xq and 
f{x) = 1 ioT X > Xq. The corresponding action functions 
are Soff{x) = x\n.{x/aQ) — x, and Son{x) = xhix — x. 
Therefore, the (logarithm of the) MSTs are given in the 
leading order of fV 1 by the accumulated action be¬ 
tween the corresponding stable metastable state and un¬ 
stable fixed point [S2] 

In Toff^on — fV[-S'(xo)-5'(ao)] = N ^xoln ^-a;o + ao^ , 

( 20 ) 

and Torinof f coincides with Tof f^on upon replacing ao by 
1. For brevity, below we only present the results for the 
off —>■ on switch. All the results related to the on —>■ off 
switch are identical upon replacing ao by 1. Note, that in 
the absence of EN, Toff^on and Ton-^f are comparable 
when xo = (1 — ao)/ln(l/ao). Eq. (20) can be simplified 
in the bifurcation limit xq — ao ^ ao. Here, we find 
In~ {N/ {2ao)){xo - ao)^ (SnUSOj- 

Now, we add negative binomial EN to the 
protein’s degradation rate. In the white-noise 
limit, integrating over the momentum 0 with 
/(x) = ao for X < xo, the action function becomes 
Soff{x) = l/(Urc) {Oao(x)-I-ln[UrcX - 1-I-Hao(x)] 
+VtcX [In ((x/2)(l — VtcX + Oa„(x))) — 1]}, where 
ilag{x) = y^(l — TeVx)'^ + AaoTcV. Therefore, the MST 
reads 

In Toff^on - N[Soffixo) - Soffiao)]. (21) 

In the bifurcation limit xo — ao ao, we find 

luToff^on ^ [iV/(2ao)](xo - ao)^/(l -b aoVre) [48]. 

In the adiabatic limit, we need to optimize the cost of 
switching from one metastable state to the other given 
noise magnitude f, against the probability of choosing 
noise magnitude f |15|. To do so we use Eq. (|T^, where 


the upper integration limit is the unstable fixed point xq, 
and for the off on switch the lower limit is the stable 
fixed point given noise magnitude p, g{f) = ao/f. As a 
result, the cost function (13) is a function of ^ only, and 
reads $o//(C) = xq [In(xo^o) - 1] + ao/C + (C - - 

1)/U. Therefore, in the leading order, the MST reads 

^ ( 22 ) 


In T, 


off- 


is the 


where = \ 1 — Uxo -b \/ (Uxq — 1)^ -b 4Uao 

saddle point satisfying = 0. 

Figure Qa-b) compares our analytical theory with 
stochastic simulations and the numerical solution of 


Hamilton equations (Bl) (see Methods section). Good 


agreement is seen in the white noise limit (similar agree¬ 
ment is seen in the adiabatic case), where the numerical 
solution of the Hamilton equations allows us to explore 
parameter ranges that are inaccessible by stochastic sim¬ 
ulation due to the long MSTs. Moreover, as can be seen 
in Figure l^c) the underlying Hamilton equations (Bl) 


capture the correct dynamics well into intermediate cor¬ 
relation time ranges, Tc = 0{1), where the white noise 
approximation breaks down. 

















Finally, we were interested in studying the effect of 
EN on population dynamics. We used Eq. (22) to cal¬ 
culate the relative fraction of the population in the on 
vs off state as a function of both the positive feedback 
threshold xq and of V, shown in Figure |^d). With zero 
or low EN the behavior of the population with regard 
to Xq is very homogeneous. Only for a very small range 
of Xq values is a macroscopically bistable population ob¬ 
served {e.g., at least 1 part in 100). If Xq is not tuned 
very precisely, no heterogeneity is observed. As the EN 
increases, though, the range of macroscopic bistability 
increases dramatically. For V > 5, well within the range 
of EN observed in biological systems, the population ex¬ 
hibits macroscopic heterogeneity across the entire range 
of a;o sampled. This effect is less pronounced, but still 
present for white EN (Figure S7, see [55]1. 



FIG. 4. (Color online) Comparison of mean switching times 
from analytical theory (lines) given by Eq. ( |21[ ), numerical 
solution of the Hamilton equations (o’s), and stochastic simu¬ 
lations (x’s) for the self-promoting gene model with N = 750, 
ao = 0.63, and xo ~ 0.80. (a) The MST from the off to 
the on state vs EN strength for white noise Tc = 0.1. (b) 
The MST from on to off for Tc = 0.1. (c) The MST vs EN 
correlation time for (Jex = 0.0365. (d) A heat map showing 
the relative probability for the system to be in the on vs off 
metastable state [using Eq. (22 l] according to the position of 
the barrier xo and the relative strength V of the EN. Here 
the EN is taken to be adiabatic. 


IV. CONCLUSIONS 

We have presented a new formalism for studying EN in 
gene expression circuits, allowing us to quantify how EN 


affects the PDFs, variances, and MSTs in various genetic 
circuits. EN is likely present in multiple forms and at 
multiple time scales [15] in a majority of processes in 
living cells. Understanding how EN alters the dynamics 
of gene networks is key for developing detailed models of 
genetic and regulatory processes. 

Our results from studying the effect of EN on simple 
genetic motifs have shown that EN has a dominant in¬ 
fluence on the system’s behavior. Analyzing experimen¬ 
tal single-cell distributions without accounting for extrin¬ 
sic noise is unlikely to provide meaningful interpretation. 
More work needs to be done to experimentally character¬ 
ize the details of the extrinsic fluctuations that cellular 
reactions rates experience. 

Finally, EN seems to provide a distinct advantage for 
populations wishing to use bistability as a bet hedging 
strategy. Without EN, the parameters needed to have 
meaningful population heterogeneity in a given condition 
are exponentially sensitive. With EN, populations are 
able to explore a variety of states across a wide range of 
parameters. 


V. METHODS 

A. Monte Carlo simulations of the auxiliary circuit 

Monte Carlo simulations with negative binomial EN 
were performed using two auxiliary species, Oi and 02 , 
to model the EN. They can be thought of in terms of a 
generic mRNA and protein, respectively. The dynamics 
of these species are given by two birth and two death 
processes: 


0 —Oi, 


ai 


0, 


Oi —> ai + 02, 

d2. 


02 


0 , 


Here, 0 is a symbol for the empty set, i.e., a species 
is created from nothing or destroyed into nothing. The 
rates are given by: 


Vl 


di 


1 K 
Tc Kal^ - 1 ’ 
ui 

7 

Tc 


V2 


d2 


- 1) 
Tc 

1 

Tc' 


K is the mean copy number of species 02 , is the 
desired EN strength, and Tc is the desired EN correla¬ 
tion time. To remind the reader, the negative binomial 
parameters, a and /3, defining the distribution (A7) are 
related to Vi and V 2 in the following manner: vi = a/xc, 
and V 2 = ojfi/Tc, where a = KjiKu^^ — 1) and /3 = 
ATcTgj, — 1. In all simulations w was set to 100. To avoid 
negative rates, one must have Ka\^ > 1. Therefore, the 
value used for K limits the lower bound of the EN that 
can be simulated using a particular set of parameters. In 
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the limit as —>■ 1 the variance approaches the Pois- 

sonian variance of a birth death process centered on K. 
Larger K allows for smaller EN to be simulated. Outside 
of this limitation, K has no influence on the EN proper¬ 
ties. It does, however, influence the computational effi¬ 
ciency of the simulation. Both larger and smaller values 
of K increase the runtime of the simulations. We used a 
value oi K = 20,000 throughout this work, which leads 
to a lower bound for the EN studied of aex > 0.007 and 
provides reasonable runtimes. 

Figure SI, see [SS], shows that the mean and variance 
of the auxiliary species are as expected during the sim¬ 
ulations. Figure S2, see ED, shows that the auxiliary 
species has the expected autocorrelation time. All sim¬ 
ulations were performed using the standard Gillespie al¬ 
gorithm m using the Lattice Microbes software 

The auxiliary species 02 was coupled to a reaction to 
be fluctuated by including it as an additional species par¬ 
ticipating in the reaction and adjusting the reaction rate 
constant such that the mean equals the original value. 
For example, to model a fluctuating birth rate for pro- 

N 

tein n with mean copy number N, the reaction 0 —> n 

becomes 02 —-—> n + a 2 . Similarly, to model a fluc- 

11 1 1 1/^ 

tuating death rate n —> 0 becomes n + a 2 - 02 ■ In 

the above equations, 02 appears on both sides to indicate 

that it is neither created nor destroyed by the reaction. 


B. Numerical solutions of the Hamilton equations 


linearize the Hamilton equations (Bl) in the vicinity of 
Xo„. This allows us to find the eigenvalues and eigenvec¬ 
tors in the vicinity of Xo„. Since the switching trajectory 
leaves the fixed point Xo„ along its unstable manifold, 
we are only considering the eigenvectors Vi and V2 that 
correspond to the two positive eigenvalues Ai and A2. 
(Note that the two other eigenvalues satisfy A 3 = —Ai 
and A4 = —A2 such that ^ ■ Xi = 0.) As a result, we 
take the initial direction to be v = Vi cos(a) -I- V 2 sin(a), 
where it is assumed that each of the eigenvectors is nor¬ 
malized to unity, and 0 < a < 27r. Finally, we search 
over all possible values of a until we find the best choice 
that satisfies the above initial and final conditions. Note 
that along Zo„(t), the initial conditions dtx(O) < 0 and 
dtPxi.0) < 0 are satisfied. 

This search over a is optimized by performing a bi¬ 
nary search. Each time an initial condition is chosen, the 
set of equations is solved numerically by using a Matlab 
numerical solver, and we compare the final condition to 
Xg. The search is terminated when we have sufficiently 
converged to the final condition. After successfully de¬ 
termining the trajectory we perform a numerical inte¬ 
gration in order to find the accumulated action. We do 
so by using the formula AS = jQ^[Px{t)x + p^{t)^]dt ~ 

T,f=ro^[x{U+i)-x{ti)] X [p^{U+i)+px{ti)]/2 + [^{U+i)- 
^(U)] X [p^{U+i) +p^(U)]/2, where to = 0 and = tf. 
This result gives us the logarithm of the mean switching 
time divided by N. An example can be seen in Figure 
S 8 in 1561. 


In this section we use the shooting method 
to find a numerical solution to the set of Hamilton equa¬ 
tions (Bl) in the case of arbitrary correlation time Tc- 
We focus on the case of the self-promoting gene, which 
gives rise to switching between metastable phenotypic 
states. Here there are three fixed points in the language 
of the deterministic rate equations, two stable points 
at Xoff = ao and Xon = 1 and one unstable point at 
Xs = xq. We are interested to numerically compute the 
trajectories, Zon{t) and Zoff(t), corresponding to the op¬ 
timal paths along which switching from the on —> off 
and off —)■ on occurs, respectively. Below we consider 
the trajectory Zo„(t), where the analysis of Zoff(t) is sim¬ 
ilar. 

Let us denote hy U = 0 and tf the initial and fi¬ 
nal simulation times, respectively. The initial condi¬ 
tion is given by Zo„(0) = Xo„ -I- (5v where Xo„ = {x = 
XomPx = 0 ,^ = l,pj = 0 ) is the corresponding fixed 
point in the 4D phase space and v is the initial direc¬ 
tion of the trajectory, see below. Here 5 is chosen to be 
small, but not too small to balance between simulation 
runtime and accuracy. The final condition is that the 
trajectory reaches the close vicinity of Xs = (xg, 0,1,0), 
namely that |zo„(</) —Xgj ^ 1. [From there, the assump¬ 
tion is that the system flows almost deterministically to 

Xo// = (oo, 0,1,0).] 

In order to find the initial direction of the trajectory we 
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Appendix A: mRNA-protein auxiliary circuit 

In this section, we derive the stationary PDF of the 
auxiliary protein. This PDF determines the extrinsic 
noise (EN) statistics of the degradation rate of the pro¬ 
tein of interest. The choice of negative binomial statistics 
used in the main text for the reaction rate seems quite 
natural. Indeed, in genetic circuits of a non-regulated 
gene, if the mRNA is short lived, the proteins’ stationary 
PDF is given by a negative binomial distribution [Mll55j . 
As a result, if this auxiliary protein affects the degrada¬ 
tion rate of our protein of interest, this rate will fluctuate 
with negative binomial statistics. 

In the auxiliary mRNA-protein circuit, mRNAs are 
transcribed at a rate a/xc and degrade with rate w/tc, 
while proteins are translated at a rate ujj3/Tc and degrade 
at a rate 1 /tc, which insures that the correlation time 
of the auxiliary proteins is We assume a short-lived 
mRNA such that w 3> 1. The master equation describing 
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the probability to find m mRNAs and k proteins satisfies: 

Pm,k — Pm^k) [(^ “f 

Tc Tc 

H- — — {Pm,k-1—Pm,k) H-[(^ + l)^m,/c+l ~ fc^rn,fc]-(Al) 

'^c '^c 

We denote the auxiliary mRNA and protein concentra- 
tions by z = vnjK and ^ = k/K, respectively, where 
K = a/3 is the auxiliary protein’s abundance. We now 
use a dissipative version of the WKB approximation, see 
e.g., Refs. [S^l [Ml IS3 ■ Employing the WKB ansatz 
Pm,k = P{z,^) ^ arrive at a Hamilton- 

Jacobi equation H = 0 with Hamiltonian 

- l)+ujz{e~P^ - 1 ) 

+/3wz(eP« - 1) + - 1), (A2) 

where Pz = dzS{z,^) and p^ = d^S{z,^) are the associ¬ 
ated mRNA and protein momenta. This yields the fol¬ 
lowing Hamilton equations 

z = ~ wze ~^‘, 

Pz = -w(e“P" - 1) - j3uj{P’^ - 1), 

^ = Pujze^^ — , 

p^ = l-e-PP (A3) 

For cli ^ 1, the mRNA lifetime is short compared to that 
of the protein. In this case, z and Pz equilibrate much 
faster than ^ and p^, and we can adiabatically eliminate 
the mRNA species [Ml | 66 ]. As a result, putting i = 
= 0 we find z = z{^,p^) and Pz = Pz(^,P^)- Plugging 
this into the Hamiltonian ( |A2| ) we arrive at the reduced 
Hamiltonian for the auxiliary protein only [21j 


Hr{^,Pi) = ^ 


- 1 


+ ae-P^-l), (A4) 


_l + /3(l-eP«) 
which effectively includes mRNA fluctuations. Solving 
the Hamilton-Jacobi equation P[ri^,P^) = 0 we find 


pj = ln[(l+/3)e/(l+/3C)], 

and thus, the action becomes 


(A5) 


5(e)=ein 


(l + /3)t 


- ^ln(l + C/3)- 


Ll + /3n 

As a result, the stationary PDF to find k copies of the 
auxiliary protein is given by P{k) ^ g--P^lSik/K)-S(i)] 
[see Eq. ([5 ) in the main text], where K = a/3 is the pro¬ 
tein abundance. This distribution, when properly nor¬ 
malized, coincides at fc ^ 1 with the negative binomial 
distribution 


Pk = 


r(a -I- k) 




1 


r(A: + l)r(a) V/S + iy V^ + iy ■ 

Interestingly, these results can give us insight on the 
mRNA fluctuations that are implicitly incorporated in 
the protein-only model [Eq. ( |A4| )[, after eliminating the 
fast mRNA variable. Indeed, by comparing Eq. (A5) 
with the momentum in the protein-only model [Eq. 0 T: 


we find that mRNA fluctuations emanating from this 
unregulated mRNA-protein circuit can be effectively ac¬ 
counted for by taking a protein-only model with a modi¬ 
fied production rate f{x) = (1 -I- (3x)/{l j3). This pro¬ 

duction rate becomes I in the limit of small burst size 
/3 —>■ 0, but becomes ^ + (I — 0//3 in the limit of large 
burst size /3 ^ 1, which yields a much wider distribution 
with variance N(3 ^ N. Importantly, this modified pro¬ 
duction rate gives rise to a protein PDF that coincides 
with the negative binomial distribution at fc ^ 1 . 


Appendix B: Analysis of the Hamiltonian combining 
IN and EN 


In this section we will derive the stationary PDF of 
the proteins of interest for generic production rate f{x), 
where the degradation rate fluctuates due to EN with 
negative binomial statistics and correlation time Tc- In 
order to do so, we will analyze the Hamilton equations 
emanating from Hamiltonian Q in the main text. In par¬ 
ticular, we will find approximate solutions for the protein 
PDF in the limits of short- and long-correlated EN. 

Using Hamiltonian Q the corresponding Hamilton 
equations read 


X = f{x)eP^ - , 

P 

p, = -nx){eP- - I) - ^{e-P- - I) 
P 

~ peP^ e~Pi^ 

r,[l +P/I-ePqY~~' 


Pi 


- 1) 
Tc 


X 

p 


{e-P- 


!)• 


(Bl) 


To remind the reader, p = K/N is the abundances 
ratio of the auxiliary protein and protein of interest 
and C = pC = is a rescaled noise variable, while 
a = — 1/K) and /3 = Ka\^ — I. Hamilton equa¬ 

tions (Bll can be solved numerically for any value of 
Tc, see Methods section. This numerical solution pro¬ 
vides the statistics of interest in the leading order, and 
is far more efficient than performing numerical Monte- 
Carlo simulations, especially for short-correlated EN, see 
main text. Importantly, the cases of fast and slow dy¬ 
namics of the EN can be studied analytically, see below. 


1. White-noise limit of EN 


In the white-noise limit, Tc <C I, the dynamics of the 
auxiliary protein is fast. As a result, ^(t) and equi¬ 
librate fast compared to x and and we can look for 
slowly-varying x and p^ dependent solutions of the third 
and fourth Hamilton equations (BI). This yields in the 
leading order of Tc I 


ffid = 1 - 2(1 - e-P-)VxTc 


(B2) 
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where we have defined V = as the ratio between 

the relative EN and IN variances. Note, that the value 
of = 0{tc) <C 1 does not enter the equations for 


X and px- As expected, Eq. (B2) as well as the re¬ 


sults below are independent of the arbitrary choice of 

the auxiliary protein abundance K. Plugging 

from Eq. (B2) into the first and second Hamilton equa¬ 


tions (B1 1 we arrive at an effective ID white-noise Hamil¬ 
tonian lainTiinH] 


H{x,px) = — l)+x{e ^'^—l)+x‘‘{e l)^Hrc. 

(B3) 

Solving the Hamilton-Jacobi equation H = 0 we find the 
momentum 


Px = In 


2f(x) I 


\ — VtcX+\/ (VtcX—I)'^ + 4:V f{x)i 


(B4) 


The PDF can be formally found by integrating Eq. (B4) 


to find the corresponding action S{x) = f Px{x')dx'\ 
and by using Eq. ([^. 


2. Adiabatic limit of EN 

In the adiabatic limit, Tc ^ 1, we can assume the EN 
is almost stationary. As a result, the stationary PDF of 
the proteins satisfies 


Pr,. = 


P{OP{nm^. 


(B5) 


where P(ji\^) is the probability to find n proteins given 
noise magnitude and P’(^) is the probability to find 
EN magnitude For simplicity we will take the EN 
to be gamma distributed, P(^) = P~°‘/T{a) 

Here a = ^l(j'lx ^nd j3 = [i/K ~ which guarantees 
that the mean is 1 and the variance is tTg,^. As can be 
checked, the gamma distribution becomes a good approx¬ 
imation of the negative binomial distribution when K is 
sufficiently large. 

With these values of a and /3, we find 


P{0^ 


I 






7T(7t 


(B 6 ) 


which holds as long as a ex <1- As a result, the PDF to 
find n proteins [Eq. (B5)] becomes 


Pr,. = 


I 


where 


'-00 C\/27rcr2 


(B7) 


P{n\^) = Ae 


-N f";,. In -M^dy 
•’gd) /(h) " 


and A = A(^) is a normalization constant. Here, we have 
used the fact that given the momentum along the op¬ 
timal path (zero-energy Hamiltonian) satisfies Px{x,^) = 
ln[x^//(a;)], and the fixed point given noise magnitude 
^ satisfies the equation f{x) = x^ and is given by 

= aiO- 


To proceed, we rewrite the integral in Eq. (B7| as 


= 


J — oo S 


(B 8 ) 


where B contains the preexponential factors including all 
normalization constants and 


4’(a;,C) = [ In 


f{y) 


dy 


V 


(B9) 


is the cost function that we need to optimize, see main 
text. Now, we use the fact that N 1 and employ 
the saddle-point approximation. The saddle point is ob¬ 
tained at = 0 , which yields the following alge¬ 

braic equation 


5$ ^ X- g{0 

- 


I 1 _ n 

4^ ■ ^ ’ 


(BIO) 


where we have used the Leibniz integral rule when dif¬ 
ferentiating Eq. (B9), and g{^) is defined above. Solv¬ 
ing the equation E[a; — 5 (^)] -f ^ — 1 = 0 for ^ yields 
the optimal noise magnitude ^*(a:). Plugging ^*(x) into 
$(x,^) we find the PDF in the adiabatic limit, which is 
given by Eq. (12) in the main text, where = 

[1 - vg'm/m- 

The variance of this PDF can be explicitly calculated. 
It is give n by the second derivative of 4)[x,^ = ^*(x)] 
[Eq. (B9) when plugging ^ = ^*(x)] with respect to x, 
evaluated at x = x» 


Nvar ^ = 


d2$[x,^ = ^*(x)] 




(Bll) 


where x* is the unperturbed fixed point (with ^ = 1 ) 
satisfying x* = /(x»). To carr y out this calculation ana¬ 
lytically we need to solve Eq. (BIOI and find the optimal 
noise magnitude We recall that the variance is cal¬ 
culated in the vicinity of the unperturbed fixed point 
X ~ X*. Let us assume a-priori that [^ — 1| <C 1 in the 
vicinity of x ~ x*. Then, we can expand g(^) in the 
vicinity of ^ = 1, g{i) ~ 5 ( 1 ) -b g'{l)il - 1). However, 
since g(l) is the solution of the equation x^ = /(x) at 
^ = 1, we have g{l) = x». Therefore, we have 


g(^) ~ Xh. -bff'(l)(g- 1). 


Plugging this into Eq. (BIO) we find 


^*(x) ~ 1 -I- 


E(x — X*) 


(B12) 


(B13) 


Vg'{l)-1' 

This verifies our assumption that 11 — (x) | <C 1 a s long 

as X is in the close vicinity of x*. Now, using Eqs. (B12) 


and (B13) in Eq. (Bll), performing the differentiation, 
and evaluating the result at x = x», we find the observed 
variance to be 


— 

^obs 


Nx,[Vg'{l)-lY 


[/'(x*)-1][2E5'(1)-1]-Ex*’ 
where we have used the fact that x* = /(x*). This ex¬ 
pression can be further simplified if we recall that g{^) 
satisfies f[g{0]/9i0 = Differentiating this with re¬ 
spect to evaluating the result at ^ = 1 , and using the 
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fact that 5 ( 1 ) = a;*, w e obt ain (/'(I) = a;*/[/'(a:*) — 1 ]. 
Plugging this into Eq. (B14) we arrive at the final result 


^2 _ 

^obs 


Nx^ 


1 - /'(a;*) 


1 + 


Pa;* 


(B15) 


1 - /'(a;*)_ 

Note, that throughout these calculations we have as¬ 
sumed that the mean of the PDF remains at a: = x*, 
and calculated the variance accordingly. This assump¬ 
tion is accurate as long as the EN magnitude is not too 
strong, cTg^, 1, see Figure S3(c) and Figure S4 in [5S] . 
which is within the range of EN observed in biological 
systems. 

Yet, for very strong EN, the mean of the PDF shifts to 
the right, due to the nonlinear dependence of the fixed 
point on the degradation rate, and due to the correspond¬ 
ing slowly-decreasing right tail of the protein PDF. We 
will now show this explicitly in the case of the unregu¬ 
lated gene, for which f{x) = 1. Let us assume = 0{1) 
such that V = = 0{N) 1. In this strong-EN 

regime, the PDF is approximately given by 

P{n) ~ C 7 \/^ev(i-w-in w)^ (b16) 

V n. 


where we have used Eqs. ( |12[ ) and (131 in the main text, 
with g{^) = 1/^ and ^*(x) ~ 1/x, and C = 

In order to calculate the mean of this PDF we use the 
equality (n) = Doing so, and using the saddle 

point approximation, we find 

{n)^N{l + 3alJ2). (B17) 

This result for the PDF mean in the case of EN in the 
degradation rate agrees well with simulations, see Figure 
S3 and Figure S4 in [5S] . 


Appendix C: The case of EN in the prodnction rate 


In this section we consider EN in the production rate 
rather than in the degradation rate. We show that while 
the resulting protein PDF in this case differs from the 
case of EN in the degradation rate, the variance of the 
PDF coincides in the two cases, in both the white- and 
adiabatic-EN limits. 

We again consider EN with a negative binomial statis¬ 
tics and correlation time Tg. Our starting point is the 2D 
Hamiltonian which encodes the stochastic dynamics of 
the protein of interest under the influence of EN. Here, 
instead of EN in the degradation rate we have EN in 
the production rate in the form f{x) —>■ ^/(x), where ^ 
satisfies (^) = 1 , and fluctuates with negative binomial 
statistics. As a result, the Hamiltonian in the case of 
EN in the degradation rate, gives way to 


P 


P 

PTc 


1 


I + /3(1 - 


- 1 


I) -|- x(e — I) 
^(e-^’5-1), (Cl) 


where ^ and p are defined above. At this point, we can 
repeat the calculations done above for EN in the degrada¬ 
tion rate. In the white noise limit we find the momentum 
to be 


Px = In 


VTgf{x) -1 + (VTcfix) - I)^ -b AVTc. 
2VTgf{x) 


(C2) 

from which the PDF can be calculated via Eq. ([^, with 
S[x) = p p{x')dx'. This PDF does not coincide with 
the case of EN in the degradation rate [compare Eq. ( |C2[ ) 
with Eq. (B4)], but for weak and moderate EN, the PDFs 
are indistinguishable, see Figure S3 (a) in (SS]. Differ¬ 
entiating the momentum with respect to x we find the 
observed variance to be 


= iV^"(x*)-^ = 


A^x*(l -I- x^Vtc) 


(C3) 


_ 1 - 

which coincides with the variance when EN is in the 
degradation rate. 

In the adiabatic case, we again need to calculate the 
integral 


Pr, 


J —00 s 


where B contains the preexponential factors including all 
normalization constants. In this case, the cost function 
$(x,^) takes the form 

y , C-ln?- 1 
>9(0 U{y) 


= [ In 

doff) 


dy 


V 


(C5) 


Note, that the only difference between this equation and 
Eq. (B9) is that here ^ is in the denominator of the In 
function, instead of the numerator. In addition, in this 
case X = g(^) solves the equation ^f(x) = x. Now , we 
use the fact that A^ 1 and solve the integral (C4) 
via the saddle-point approximation. The saddle point 
is obtained at 9j$(x,^) = 0 , which yields the following 
algebraic equation 


giO-x 


I 

V 


I 


= 0 . 


(C 6 ) 


^ V 

Solving the equation H[(7(0 ~ 2 ;] + C ~ 1 = 0 for ^ we 
find ^*(x), which allows finding the PDF according to 
Eq. (12), see Figure S3(b) in [56]. This figure emphasizes 


the lack of coincidence between the PDFs in the cases of 
adiabatic EN in the production and degradation rates. 

The variance of this PDF is given by Eq. (Bill, 
where x = x* is the unperturbed fixed point satisfying 
a;* = /(x*). To carry out this calculation analytically we 
need to solve Eq. ( |C6| ) and find the optimal noise magni¬ 
tude We recall that the variance is calculated in the 
vicinity of the fixed point x ~ x*. Assuming a-priori that 
1^ — 1| ^ 1 in the vicinity of x ~ x*, we take Eq. (C6) 
and expand g(^) to first order in ^ around ^ = 1. By 
doing so, and using the fact that 5 ( 1 ) = x*, we have 
g(0 - 1), which yields 


^*(x) ~ 1 -I- 


V(x — X*) 
Vg'il) + 1' 


(C7) 
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Indeed, we find that |1 — ^*(a;)| <C 1 as long as x is 
in t he cl ose vicinity of x*. Now, we p lug <I>(a:,^) from 
Eq. (C5) and ^ = C*{x) from Eq. (C7) into Eq. ( |B11 ). 
Performing the differentiation twice with respect to x, 
plugging X = X* = /(x*), using the fact that g{^) satisfies 
/b(?)]/ff (0 = 1 /C) and evaluating the result at ^ = 1 
which yields ^'(l) = x*/[l —/'(x*)], we find the observed 
variance to be 


^obs 


Nxt 


1 + 


Ex* 


1 - [ 1 - 


(C 8 ) 


This result coincides with the variance in the case of EN 
in the degradation rate, see Figure S3(d) in [5S] . 


Appendix D: The case of Ornstein-Uhlenbeck EN 

In this section we consider EN with different statistics. 
We take Ornstein-Uhlenbeck (OU) extrinsic noise with 
mean (^(t)) = 1 and variance (C(^)C(^O) = 
with correlation time Tc > 0. Note that in our previous 
work [35] on the self-regulating-gene model, we have al¬ 
ready used the OU noise when modeling EN. Here we 
develop a different and more generic formalism allowing 
to go beyond the bifurcation limit done previously, and 
to treat EN of arbitrary strength. Notably, EN with such 
statistics can give rise to zero or even negative reaction 
rates for sufficiently strong EN, which can cause, e.g., the 
divergence of the mean |69j . As a result, in our derivation 
below we implicitly assume that the noise statistics has 
a cutoff such that the reaction rates are always positive 
real numbers. 

The OU process satisfies the following Langevin equa¬ 
tion 


C = - (C - + V 2 cr 2 ^/rc 77 (t), (Dl) 


dynamics [Eq. ([^] . This yields a 2D master equation for 
the probability P{n,k,t) to find protein copy number n 
and noise copy number k, at time t. Similarly as in the 
case of negative binomial EN, using the WKB ansatz for 
the stationary PDF, Pn^k ^-NS{n/N,k/N) ^ arrive at 
a Hamilton-Jacobi equation H = 0 with a Hamiltonian 


H^'^^'>{x,Pa:,i,p^) = - 1) + —(e 

^ P 

(2p^U - f + p) ^ (2p^U + |-p) ^^ 


2Tr. 


2Tc 


- 1) (D3) 

- 1 ), 


where as before V = ^ and p = K/N, while 

Px = dxS, and p^ = d^S are the associated momenta. 
This Hamiltonian encodes the stochastic dynamics of the 
protein of interest when its degradation rate fluctuates 
with OU noise. 

In order to proceed, we can write down the correspond¬ 
ing Hamilton equations 


X = f{x)e^^ - 

P 

Px = - 1) - - I) 

P 

X pP{ _ p-Pl 

c = - e + p) - + e - p) 

pr = ^{en - e-P£) - -{e-P- - I). (D 4 ) 

« 2tc p 


Similarly as in the negative binomial case, in the white- 
noise limit, Tc 1 , ^(t) and p^{t) equilibrate fast com¬ 
pared to X and Px- As a result, we can look for slowly- 
varying X and Px depende nt s olutions of the third and 
fourth Hamilton equations (D4), which yields in the lead¬ 
ing order of Tc 


where p(t) is white noise {g{t)r]{t')) = S{t — t'). Here 
r]{t) can be defined as the dt —>■ 0 limit of the temporally 
uncorrelated normal random variable with mean 0 and 
variance 1/dt. The stationary statistics of this noise is 

In order to go beyond the bifurcation limit |48| , we are 
interested to describe the OU process via a discrete birth 
death process describable by a master equation. Defining 
k = as the noise “copy number” in the OU process 
where AT ^ 1 is an arbitrary large number, the master 
equation describing the probability Pk to find EN copy 
number k satisfies: 


Pk = ^k-lPk-l + Vk+lPk+1 — (A/c + Vk)Pk- (D2) 
Here \k = l/{2Tc){2.K‘^a‘lx — k -\- K) and Vk = 

l/(2Tc)(2A'^crg,j, + k — K) are the birth and death rates, 
respectively. One can check that using these birth and 
death rates one recovers the Langevin equation for the 
EN “copy number”: k = —{k—K)lTcP\/2JPPaPjVc g{t), 
which corresponds to Eq. (Dl I with ^ = k/K. 

To study the interplay between IN and EN, we com¬ 
bine the EN dynamics [Eq. (D2)] with the underlying IN 


= I - 2{1 - e-P’^)VxTc. 


(D5) 


Note, that the value of p^~^^ = 0 {tc) ^ 1 does not enter 
the equations for x and Px- Also, one can see that this 
result as well as the results below are independent of 
the arbitrary choice of K. Pluggin g thi s effective noise 
into the first of Hamilton equations (D4) we arrive at an 
effective ID white-noise Hamiltonian 


H{x,px) = f{x){eP^— l)+x{e —l)-|-x^(e — 1 )^Utc, 

(D 6 ) 


which coincides with Eq. (B3|. As a result, in the white- 


noise limit the protein PDF under OU extrinsic noise 
coincides with the case of negative binomial EN. In par¬ 
ticular, the variance in this case coincides with Eq. ( |To| . 
More generally, this indicates that in the white-noise 
limit, the choice of EN statistics does not affect the PDF 
in the leading order of Tc <C 1 . 

In the adiabatic regime, 1, similarly as 

in the negative binomial case, we can use Eq. (B5) 
with P(^) = Xj^2'Ka\x: and P{n\^) = 

where x = g{^) solves the equation 
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f{x) = x^. As a result, we arrive at Eq. (B8), where 
here, the cost function satisfies 


= [ In 

Jgii) 


f{y) 


dy - 


2 E 


(D7) 


Now , we use the fact that N ^ 1 and solve the integral 
in Eq. (B8) via the saddle-point approximation. The sad¬ 


dle point is obtained at d^^{x,^) = 0. Using Eq. (D7), 
this yields the following algebraic equation 


x-g{^) f - 1 


= 0 . 


(D8) 


^ ' V 

Solving the equation U[a; —g(^)]-|-^(^ —1) = 0 for ^ yields 
the optim al noise magnitude Using this result and 

Eq. (D7), we find the PDF according to Eq. (12) in the 


main text. Note, that the resulting PDF here differs 
from the negative binomial case, since the cost function 
here [Eq. ( |D7[ )] differs from that in the case of negative 
binomial EN [Eq. (B9)]. 


The variance of this PDF can be explicitly calculated 
by using Eq. (Bll I, with x* being the unperturbed fixed 
point X* = /(x*). Since the variance is calculated in the 
close vicinity of the fixed point x ~ x», similarly as for 
the negative binomial EN, we find the saddle point to be 


^*(x) ~ 1 -I- 


Vg'{l) - 1 ^ 


(D9) 


which coincides with Eq. (B13). Repeating the calcu¬ 


lations in the same manner as in the case of negative 
binomial EN, we find 


' obs 


Nx, 

1 - /'(x*) 


1-h 


Ux* 


(DIO) 


l-/'(x*)J' 

This result again coincides with the variance in the neg¬ 
ative binomial case [Eq. (B15)]. This indicates that to 
determine the variance of the protein PDF under EN (in 
both the white- and adiabatic-noise limits), the complete 
statistics of the EN is less relevant. The only relevant pa¬ 
rameter here is the width of the EN distribution, or its 
magnitude, given by the parameter a^x- 


Appendix E: Correction of analytical variance using 
the numerical mean 

In cases where the EN magnitude is large, the mean of 
the distribution can shift, as discussed above. In these 
cases we apply a correction to the observed variance to 
account for the change in the IN. For example, = 
N{1 + Vtc) is corrected to 

'^L = ^«(l + ^rr,), (El) 

where gobs is the mean observed from numerical simula¬ 
tions, and V = Na^x- 
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